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Abstract 


Despite the variety of methods available to solve nonlinear optimal con- 
trol problems, numerical methods are still evolving to solve these problems. 
This paper deals with the numerical solution of nonlinear optimal control 
affine problems by the interpolated variational iteration method, which was 
introduced in 2016 to improve the variational iteration method. For this 
purpose, the optimality conditions are first derived as a two-point bound- 
ary value problem and then converted to an initial value problem with 
the unknown initial values for costates. The speed and convergence of the 
method are compared with the existing methods in the form of three ex- 
amples, and the initial values of the costates are obtained by an efficient 
technique in each iteration. 
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1 Introduction 


One of the important problems in the field of engineering [13, 6], robotics 
[14, 11], biology [1, 16], and so on, is finding the best control strategy to 
use available resources optimally. In practice, most of these problems are 
formulated nonlinearly, which in general do not have analytical solutions, and 
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one must look for reliable numerical solutions for them. Generally, consider 
the following nonlinear control system: 


&(t) = A(x(t))+ Bidju(t), — t E [to ty], 
r(to) =X, «z(t>) = zz, 


(1) 


where x(t) and u(t) are n- and m-vector of state and control functions, re- 
spectively. The aim is to minimize the objective functional 


Java) == fa QU a(t) + uP )R()ule))at (2) 


to 


subject to the dynamical system (1). Moreover, Q(t) is a positive semi- 
definite real matrix, and R(t) is a positive definite real matrix. 


There are two approaches to solving this problem. One is the direct 
method of parameterizing the state and control functions and solving a non- 
linear optimization problem. The other is the indirect method of using opti- 
mality conditions, derived by the Pontryagin maximum principle (PMP) or 
Hamilton—Jacobi-Bellman equation [3, 15]. The main topic of our discussion 
in this paper is the optimality conditions derived by PMP, which results in 
the following two-point boundary value problem (TPBVP) as 


&(t) = A(a(t)) + Bult), 
At) = —Q(ta(t) — AT (a(t) A(t), 
x(to) = Xo, u(tf) = xe, (3) 


where (t) is an n-vector of co-state functions, AZ (x(t)) = VzA7 (a(t)), and 
the optimal control is given by 


u(t) = —R-! BT (t)X(t). (4) 


In 2012, Berkani, Manseur, and Maidi [4] used the variational iteration 
method (VIM) to solve this problem, for which the control variable can be 
calculated in terms of state variables. They then turned the problem into a 
variational problem and solved it with VIM. It is noteworthy that the VIM is 
an analytical-approximate method for solving initial value problems (IVPs), 
which was introduced by He [8] in 1999, and many researchers developed 
the method to increase its efficiency and validity in recent years [20, 5]. For 
solving (3), the author and Effati [24] used a combination of the shooting 
method with VIM. However, the time-consuming and complex integration of 
nonlinear functions reduces its efficiency, at each iteration. 


Matinfar and Saeidy [17] applied a modified VIM (MVIM) that uses 
He polynomials. Simultaneously, Saberi Nik, Effati, and Yildirim [22] used 
a differential transformation method along with differential transformation 
polynomials to solve (3). A comparison between the existing methods, 
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MVIM, Homotopy perturbation method (HPM), and Adomian decomposi- 
tion method was done in 2016 by Jafari. Ghasempour, and Baleanu [9]. 

Mirhosseini-Alizamini and Effati [19] extended the VIM for optimal delay 
control problems. In 2020, the author and Effati [25] used the successive 
approximation method (SAM) to solve (3), which has a very good computa- 
tional speed but does not guarantee the convergence for problems with large 
domains. 

An improvement of VIM, called Interpolated VIM (IVIM) introduced by 
Salkuyeh and Tavakoli [23], interpolates the VIM with the well-known spline 
functions, which make the method suitable for approximating all types of 
L?{to,t] solutions, especially for control and state functions. This property 
motivates us to apply the IVIM for solving optimal control problems. This 
method has low computational time in comparison with VIM since it uses an 
approximation for integration. It also maintains accuracy and convergence. 

This paper is organized as follows. In Section 2, the basic ideas of IVIM 
are briefly described. Then, the application of IVIM for solving nonlinear op- 
timal control problems (NOCPs) is explained in Section 3. Several numerical 
examples are illustrated in Section 4, which show the accuracy and efficiency 
of the proposed method compared to the existing ones. 


2 Basic ideas of the interpolated variational iteration 
method 


As stated earlier, the VIM is an analytical-approximate iterative method 
for solving IVPs [8]. To explain this method, consider the following one- 
dimensional IVP 


fe eae teE (0, 7], (5) 
where y(-) € R and f satisfies the Lipschitz condition: 
f(t, y(t) — FE IO) < Ly) — 9@)I, (6) 


for any t € [0,T] and for some Lipschitz constant L > 0. We make the 
following iterative formula, which is known as the VIM main formula: 


te (d+ f Mette Hawa. a0 


where A(s,t) is known as a Lagrange multiplier [2, 28], and yo(t) should be 
chosen such that it satisfies yo(0) = y?. The VIM formula produces a se- 
quence of functions {Yym(t)}, which converges to the exact solution of (5), 
under suitable conditions and does not require any linearization or pertur- 
bation. Due to the attractiveness of the method for calculating accurate 
solutions, many researchers have used it to solve difficult problems such as 
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dispersive optimal solitons [7], coupled Burgers’ equation [26], fractional op- 
tion pricing [18], fractal Fisher’s equation [27]. However, the integration of 
nonlinear function f(t, y(t)) in all these researches may greatly reduce the 
speed of the method or lead to the failure of the method. To overcome this 
problem, first, we apply the integration by parts for (7) to get: 


imal) = Gnl)= | Hn(s, t)ds, (8) 
where 
Gin(t) = (1+ A(E,#))¥m(t) —-A(0, #) 4m (0), 
and 


Hy(5,t) = (6, 8)tin(s) + Ms, 8)F(6,4m(8)) 


Now, we interpolate y(t) by first order B-spline functions [12], at equidis- 


tance points t;, 7 = 1,...,n on [0,7]. For any i = 2,3,...,n —1, let us 
define a ; ‘ 

— ti-1 it1 — 

B;(t) a Micato) + X[titepr) (t), 
h h 
and ty ag 
~ bn—1 
B,(t) = | Xitn=aitn) M4); 


where y.4 is the characteristic function defined on a set A. Its value is 1 on 
A and is 0 otherwise. It is obvious from the definition that 6;(t;) = 1 and 
B,(t;) =0, for 7 # i. Then, we employ the interpolant 


Ym(t) = 2 Ym (ti)Bi(t). 


Thus, §m(ti) = Ym(ti), for i = 2,...,n. Now, we substitute the interpolant 
Gm into (7) and after some simplifications, we arrive at 


Ym-+1(ti) = Gm-+1 (ti) 


This is the main formula of IVIM, which overcomes the nonlinear integration 
of VIM and will be used as a tool for solving NOCPs in the next section. 

Now, we state the convergence of VIM and IVIM sequences as the next 
two theorems, which were proved in [23]. 


Theorem 1. Let f be a continuous function on [0,7]. If ym(t) € C*((0, 7) 
for m = 0,1,2,..., then the VIM sequence (7) converges to the exact solution 
of (5). 
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Theorem 2. Let f be a two times continuously differentiable function. If 
the VIM sequence (7) is to the exact solution of (5), then so is the IVIM 
sequence (9). 


3 The application of IVIM for solving NOCPs 


For the sake of simplicity, we rewrite the TPBVP in (3) as 


t= rats x, A), 
\=G(t, x, d), 
a(to)= xo, a(t) = xf, (10) 


where 


F(t,x,) = A(a(t)) — B(t)R-' BT (t)X(t), 
G(t, x, ) = —Q(t)a(t) — Az (a(t))A(). 


This nonlinear TPBVP has no analytical solution in general. Thus, re- 
searchers are looking for reliable and almost accurate solutions to it. The 
technique we are studying here is IVIM, which was described in Section 2. 
Given that the above problem is TPBVP, we first convert it to an IVP as 
follows [24]: 

& = F(t,2,), 

AN=G(t,2,r),  t € [to, ty] (11) 

t(to)= 20, A(to) =a, 


where q@ is an n-vector unknown variable. This variable could be identified 
after sufficient iterations of IVIM. 


As Section 2, we create the VIM correction functional as follows: 


tmea(t) =am(t)+ f Asti) Ge. Flsante inlaid G2) 


Anea(t) = Alt) + fh Aattes) {Aon(s) — G(s. arm (s)sAna(s))} a (13) 


with xo(t) = ao, Ao(t) = a and m > 0. To determine general multipliers 
Ai(t,s) and A9(t,s), we restrict the variations for F and G, that is, 5F = 
6G =0. Therefore, the optimal multipliers become Ax(t, 8) = Ao(t,s) = —1. 
Compared with (8), consider Gi(t) and H(t) corresponding to x(t) and 


Gi m/(t) = (1 + Ai(t, t))am(t) = Ay (0, t)xo, 


Hi m(s,t) = POs, t)tm (8) + Ai(s,t)F(s8,@m(s),Am(s)). 
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Similarly, consider G2(t) and H2(t) corresponding to A(t) and 


Gom(t) = (1+ Ag(t, t))am(t) — A2(0, tha, 


OA 
Ham (3st) = “pe (5 t)#m(s) + Ao(s,t)G(s,%m(s8), Am(s))- 
For a given n € N, assume ¢; = to + (i— 1)h, i = 2,3,..., where h = “1-0 
Then, the B-spline interpolants of %m4+1(t) and Am+1(t) become 
h i 
itt) = », &m4+1(ti)Bi(t), 
i=2 
MAO) = So Amt BO), 
i=2 
where the coefficients are calculated as 
Lm+i1(ti) = Fm41(ti) 
i-1 h 
= Gi m(ti) — 3 Hy m(tr, ti) — 3 (Aim (ti, ti) + Him (ti, ti)) , 
(14) 


Am+1 (ti) = Am+i (ti) 


i-l1 
h 
= Gom(ti) —h S Ho m(tr, ti) — 3 (Hom (ti, ti) + Ho m(ti,ti)) , 
r=2 


(15) 


for any 7 = 2,3,...,n. We call these two equations as the IVIM main formu- 
las. The convergence of the IVIM is a straightforward conclusion of Theorems 
1-2, which we state in the next theorem. 


Theorem 3. Let F and G be two times continuously differentiable functions 
with respect to their first variable and let they satisfy the Lipschitz condition. 
If &m(t), Am(t) € C*([to,t¢]), then the IVIM of (14)-(15) is convergent to 
the exact solutions of (11), say x*(t,a@) and A*(t,a). Accordingly, x*(t,a*), 
A*(t,a*) are the exact solution of (10) when a* is the real root of x*(tr,a) = 
xf. 

Proof. Since F and G are continuous on [to, ty], and 2m (t), An (t) € C*([to, tf]), 
Theorem 3 guarantees the convergence of the VIM formulas (12)-(13) to 
the exact solution of (11). Besides, since F and G are two times contin- 
uously differentiable functions, Theorem 4 ensures that the IVIM is con- 
vergent to the exact solution of (11), say a*(t,a@) and A*(t,a). Clearly, 
for every a, x*(tp,@) = x and A*(tg,a~) = a. Then, from (10), we have 
x(tr) = x, = x*(ty,a). Therefore, to get the exact solution of (10), one 
should solve x*(t+,a@) = xy for some appropriate real a. 
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Theorem 4. Under the assumptions of Theorem 3, the sequences {uy,(t)} 
and {J,,} defined by 


n-1 

(L(xo, tm(tr)) + L(am(ts),0)) +h D> L(am(ti), Um(ti)), — (17) 
i=2 

converge to the optimal control and optimal objective value, respectively, 

where L(x, u) = $(a7 Qa + u’ Ru). 

Proof. Taking limit from (16) and putting a = a* yield 

* _—_ Tt — _p-lpT : 
u'(t) = im Um(t) = —R-*B’ (t) im Am(t) 


= -RBT (ON (0°), 


where u*(t) is the optimal control function, since «*(t,a*) and A*(t,a*) are 
the exact solutions of extreme conditions in (10). 


Now, taking limit from (17), as n,m — oo, and using the trapezoidal 
integration rule result in 


lim lim Jmn 
m—+oo N00 . 


= lim lim * (£(20;Um(tr)) + £2 m(tz),0)) +S Ll m( ti) (ti) 


m—>co NCO 


I 
. 
5 
pa 
8 
: 
_ 
iG 
S 
= 
ae, 
oS 
= 
ot 


= : i (Jim 2n()QW) im tm(t) + lim u;, (R(t) Jim un(t)) dt 
7 5 i. (a*(t,a*)Q(t)2*(t,a*) + u*(t)" Ru*(t)) dt 
=" 


It follows that J* is the optimal objective value and the proof is completed. 


Now, we present an algorithm to summarize our proposed method. We 
choose the stopping criterion 


<€, (18) 
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which maintains the accuracy of the objective value for any sufficiently small 
tolerance limit € > 0. 


Algorithm of IVIM for solving nonlinear OCPs 


Step 1. Let m:=1, xo(t) := 2%, Ao(t) := a and e > 0. 

Step 2. Calculate state and costate functions, 2m+1(t,@) and Am+i(t, a), 
through the IVIM formulas (14)-(15). 

Step 3. Solve em+41(ts,a) = xf for a, to get a®. 

Step 4. Find the suboptimal control u.+1(t,a*) and the suboptimal objec- 
tive value Jm+i, by (16)—(17). 

Step 5. If Em4i <€, then goto Step 6; otherwise let m := m+ 1 and goto 
Step 2. 

Step 6. Stop the Algorithm. u,+1(t) and Jm41 are the desired optimal 
control and objective value. 


Remark 1. If the state function has free final condition, then the costate 
final condition A(t) = 0 must be applied. Therefore, to calculate a in Step 
3, the equation »,,+1(ty,a@) = 0 must be used instead of tm41(t,a@) = xf. 


4 Illustrative examples 


In this section, we give two nonlinear examples to illustrate the reliability 
and efficiency of our proposed method. In the first example, we describe the 
algorithm process in detail. Also, we present the second the third example 
to compare the method with some other existing methods. Since there are 
no exact solutions for these three examples, we use criterion (18) to stop 
the algorithm and check the accuracy of the method. Besides, we use the 
Maple “dsolve” package as a reliable solution for comparison. The codes 
were developed using symbolic software Maple 18 on a machine with Intel 
Core 2 Due Processor 2.53 GHz and 4 GB RAM. 


Example 1. Consider the following nonlinear control system [24, 25]: 


& = =x" (t)sin z(t) + u(t), t € [0,1], 
x(0)=0, a(1)=0.5, 


under which the following objective functional should be minimized: 


t= [ u?(t)dt. 


Conforming to (3), we have 
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a= sot sin x(t) — 5Nt), (19) 
Mt) = —A(t)x(t) sin x(t) — 5 NE)2°(6) cosx(t), t€ [0,1] (20) 
x(0)=0, «(1)=0.5, (21) 
and the control function is u(t) = —5A(t). We substitute the final state 


condition z(1) = 0.5, by A(0) = a to get 


where a is an unknown real variable. In view of (11), F and G are the right- 
hand sides of (19)—(20), respectively. They are analytic functions with respect 
to x and A, and they satisfy the Lipschitz condition (6), with L¢ = 9/16 and 
Lg = 5/16, whenever |z| < 0.5 and A < 1. Hence, Theorem 3 guarantees the 
convergence of both VIM and IVIM formulas. The VIM formulas (12)—(13) 
become 


Lm+ti(t) =Lm(t) — | {im(s) — stl) SiN Up,(s) + 5 Am(s) ds, 


Am+1(t) =Am(t) -| {Am(s) + Am(S)¢m(s) sin m(s) 


+ yAm(s)22,(s) COS L,(s)}ds, 
xo(t) =0, Ao(t) =a, m>0. 


Now, interpolating VIM at t; = 0, 4,1, the IVIM formulas (14)—(15) give 
the following results, after the first iteration: 


x (t;) = 0, 4° jag Sin( Q), gae Sint a) a 39% Sin( Q), 
1 1 
Hs 26 3 
Ai(t;) =a, a 16% sin(7 a) 138% cos(—a), 
a Penta : a er) a sin(a) : a’ cos(=a). 
64 32 
The final state condition, x1(1) = 0.5, results in 
es a: Cn re 
; = 0. 22 
64% sin(7@) 5% 35% sin(5a) 0.5, (22) 
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Table 1: Calculation of a for n = 3 and m = 1, Example 1. 


a | Objective functional 
—0.9659445089 0.2146370242 
—1.023796574 0.2387946125 


which obtains two real roots. These real roots along with their corresponding 
objective values are summarized in Table 1. We suggest choosing the parame- 
ter a, which gives a better objective value, that is, a = —0.9659445089. This 
parameter could be a good starting point for root-finding iterative methods, 
such as Newton, at the next iteration of the IVIM. 

Table 2 shows the results of [VIM procedure up to m = 5 iterations with 
n = 10 nodes. 


Table 2: Simulation results for n = 10, Example 1. 


5 Em CPU time a 
0.22264786 - 0.000 -0.9717459824 
0.23465678 | 5.1177E-02 0.016 -0.9985940802 
0.23541609 | 3.2254E-03 0.078 -0.9997821239 
0.23513137 | 1.2109E-03 0.359 -0.9991493164 
0.23512939 | 8.4519E-06 2.266 -0.9991493164 


ap wnm eS 


To compare the IVIM with two recent methods, VIM [24] and Modal 
series [10], we summarize their results in Table 3 after five iterations with the 
same stopping criterion (18), with e = 2 x 10~°. Clearly, the IVIM converges 
rapidly, in contrast with VIM and Modal series, since it reaches a better 
tolerance limit. Also, the suboptimal control and state of the proposed VIM 
and Modal series methods are demonstrated in Figure 1. 


Table 3: Comparing results of [VIM, VIM, Modal series, Example 1. 


Method | ies is 
IVIM (m=5) 0.2351 | 8.5E-06 
VIM (m=5) 0.2353 | 1.7E-03 


Modal series method (m=5) | 0.2346 | 1.3E-03 


Example 2. Consider the nonlinear control system described by 
Ly =XQ+ 2122, 
fg = 2, +%2 + U3 +, (23) 
r1(0) = —0.8, r2(0) = 0, 


and the objective functional 
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0.54 


0.14 


t t 


Exact O IVIM 


Exact O IVIM 


Figure 1: Suboptimal state and control functions, Example 1. 


1 1 
I= 5/ (aj + 25 + u?)dt. (24) 


One can obtain the optimality conditions as 


AL = —(%1 + A1x2 — 2), 

2 = —(a2 + Ai(1 + £1) + A2(1 + 2x2), 
Ly =09 2122, 

L2 =—-2,+29+2%—2o, 

x1(0) = —0.8, x2(0) = 0, 


and the optimal control law as u = —A 9. Since x; (t) and a(t) are free at 
t = 1, we use the Remark 1 and put A;(1) = Ag(1) = 0. 


Now, due to the polynomial nature of the system and their validity in the 
assumptions of Theorem 3, we apply the proposed algorithm for « = 1073 
to get an accurate optimal objective value. It is worth noting that for some 
iterations, m = 2, solving Aim(1,a@) = 0 and Az m(1,a) = 0, results in 
different values of a. In these cases, as shown in Table 4, we clearly select 
the vector a that yields the minimum objective functional, that is, ay = 
—1.227402971 and ag = 0.5053886518. 


Table 4: Calculation of a for m = 2, Example 2. 


Qa, | a2 | Objective functional 
-1.227402971 | 0.5053886518 0.4113400952 
4.538547207 | -1.207206466 2.664175385 
1.875104138 | -1.907110672 4.484594474 
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Table 5: Simulation results for different iterations, Example 2. 


VIM SAM IVIM 
Em Time Jim Em Time Jey Em Time 
- 0.031 0.5357 - 0.015 0.8560 - 
0.4108 | 1.1E-0 | 0.031 0.7575 | 2.9E-1 | 0.031 0.4113 | 1.1E-0 | 0.015 
0.4325 | 5.0E-2 | 0.046 0.9930 | 2.4E-1 | 0.063 0.4340 | 5.2E-2 | 0.078 
1.6E-2 | 0.125 0.6750 | 4.7E-1 | 0.046 0.4413 | 1.7E-2 | 0.124 
0.4422 | 6.1E-3 | 1.544 0.3503 | 9.3E-1 | 0.094 0.4443 | 6.9E-3 | 0.296 
0.4430 | 1.8E-3 | 45.693 0.4245 | 1.7E-1 | 0.078 0.4453 | 2.3E-3 | 1.279 

- - >1800 0.6025 | 3.0E-1 | 0.078 0.4456 | 6.1E-4 | 7.504 


3 


— 
= 
+ 
A 

a 
YQ 

3 


Nooakwne gs 
So 
is 
oo 
oS 
ot 


Then, one can use these a; and az as initial guesses for solving A1,3(1, a) = 
0 and 23(1,a) = 0. This technique greatly improves system speed and 
reduces computational complexity. 


The objective values, relative errors and the CPU time of the algorithm 
iterations are summarized in Table 5. Moreover, to clarify the reliability 
of the algorithm, we compare the results with VIM [24] and SAM [25]. As 
shown in Table 5, SAM has a divergent behavior, and although it has less 
computational time than the other two methods, it does not reach the desired 
tolerance limit. Also, the VIM requires more than 1800 seconds of computa- 
tional time for m = 7 and does not reach the desired tolerance. This is while 
the IVIM reaches the desired tolerance limit after m = 7 iterations, and it 
only takes 7.504 seconds of CPU time. 

From the beginning of the algorithm, if we considered e = 5 x 107°, both 
IVIM and VIM would get the desired solutions after m = 6 iterations, but 
IVIM in 1.279 seconds and VIM in 45.693 seconds of CPU time, which clearly 
shows the speed of convergence of IVIM. 

Figure 2 shows the solutions after m = 7 algorithm iterations, in compar- 
ison with “dsolve” package of Maple as the exact solutions. 


Example 3. Consider the optimal maneuvers of a rigid asymmetric space- 
craft given by [21] 


in(t) = 2? na(t)aalt) + Palo, 
a(t) = A nWirolt) + Dual, (25) 
a(t) = - 2? 2, (t)aa(t) + F-wslt), 


with the boundary conditions 


x1(0) =0.01r/s, «x2(0) =0.005r/s, 2x3(0) = 0.001 1r/s, 
x1(100) = x2(100) = x3(100) = Or/s, (26) 
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0.74 0.8 
-0.75 ie 
0.6 

-0.76 
0.5 
= -0.77 04 
0.78 ae 
0.2 

0.79 
0.1 
-0.806 0 

0 o2 O04 06 O08 1 0 02 O04 06 08 1 


Exact O IVIM 


0 0.2 0.4 0.6 0.8 1 


Exact O IVIM 


Figure 2: The optimal trajectories and control of Example 2. 


where £1, £2, and x3 are angular velocities of the spacecraft, ui, uw2, and us 
are control torques, J, = 86.24, Iz = 85.07, and Iz = 113.59 kg m? are the 
spacecraft principle inertia. The controls u;(t),i = 1,2,3, should be chosen 
to minimize the following quadratic performance index: 


100 
Jle,1] = 5 ‘) (u(t) + ud(t) + w2(t))ae, (27) 


oe 


Also, B(t) = diag( ), Q(t) = O3x3, and R(t) = Isx3. Then, the 


eee 
optimality conditions become 
; Iz — 1. Ai(t 
#1(t) = — 2 ? x2(t)x3(t) — u(t) 
I, Ty 
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ia(t) = — 9 au(taatt) - “By, 
is(t) = — 25 an(taatt) - “By, 
A(t) = + salt) ralt) + 2 aalt)rale), 
Ja(t) = 2 a(t) u(t) + 2 cr(rald), 
Ja(t) = 2 aaltyant) + 2 ar (t)rale), 


with the boundary conditions (26), and the optimal control laws are given 
by 


, #=1,2,3, t€ [0,100]. (28) 
By using the proposed algorithm with « = 10~°, we reach the desired solution 
after m = 6 iterations. 


Table 6: Simulation results for n = 10, Example 3. 


m dni Em ay a2 a3 

1 | 0.00468914 - 0.7314359 | 0.3851119 | 0.1356719 
2 | 0.00468570 | 7.3285E-04 | 0.7453951 | 0.3581385 | 0.1272553 
3 | 0.00468711 | 2.9978E-04 | 0.7437911 | 0.3617697 | 0.1281838 
4 | 0.00468791 | 1.7023E-04 | 0.7437022 | 0.3619104 | 0.1291345 
5 | 0.00468782 | 1.9383E-05 | 0.7437007 | 0.3619099 | 0.1290549 
6 | 0.00468781 | 7.9781E-07 | 0.7437039 | 0.3619034 | 0.1290513 


The simulation results along with the estimation of unknown parameters 
a;,1 = 1,2,3, are given in Table 6. 


Table 7: The maximum error of x; (t) in comparison to HPM [21], Example 3. 


m HPM IVIM 
(Itr.) | Max Error CPU time Max Error CPU time 
2 7.3102E-05 7.328 9.7727E-06 0.031 
3 9.2612E-06 12.542 5.7522E-07 0.109 
a 1.9301E-06 18.664 4.4215E-08 0.188 
5 2.2560E-07 36.729 3.7427E-08 0.655 
6 3.1420E-08 46.401 1.6156E-08 3.994 


To compare the solutions with HPM [21], we have summarized the maxi- 
mum error of z(t) in both methods compared to “dsolve” package of Maple 
as the exact solution, In Table 7. As is evident, the IVIM is faster and 
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more accurate than HPM, so as to reach the error of 5 x 107°, HPM needs 
m = 6 iterations with 46.401 seconds of calculation time, but IVIM, needs 4 
iterations with 0.188 seconds of calculation time. 


0.010 & 0.005 & 
0.008 5 0.0044 

0.006 5 0.003 

Ye 
0.0044 0.0027 
0.002 4 0.0014 
0+ r r r r Q 04 r r r r 
0 20 40 60 80 100 0 20 40 60 80 100 


Exact O IVIM 


Exact O IVIM 


0.0010 & 


0.0008 4 


0.00024 


0 20 40 60 80 100 
t 
Exact O IVIM 


Figure 3: The optimal states x1(t),ro(t), and x3(t), Example 3. 


Figures 3-4 show the optimal states and controls of [VIM for m = 6 
iterations compared to “dsolve” solutions. The maximum errors of our pro- 
posed method are 1.6156 x 107°, 1.0231 x 107°, and 1.39963 x 1078, for state 
functions x1, 22, and x3, respectively, and 3.4660 x 10-7, 6.8377 x 10~7, and 
2.1526 x 10~", for control functions u;, u2, and u3, respectively. Meanwhile, 
only 3.994 seconds are needed to perform this operation. 


5 Conclusions 


In this paper, the interpolated VIM was used to solve NOCPs. This method 
used basic spline functions as interpolants, which is desirable for the approxi- 
mate calculation of state and control functions. This method did not require 
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-0.00856+ 
-0.008574 


-0.00858 


=~ -0.008594 
-0.008604 


-0.008614 


-0.008624 


0 2 2440 «60 80 100 


Exact O IVIM 


0 20 40 60 80 100 
t 


Exact O IVIM 


Figure 4: The optimal controls uj (¢),w2(t), and u3(t), Example 3. 


complex analytical integrations and therefore was faster than VIM. Further- 
more, the convergence range of the [VIM was wider, compared to the new 
SAM. In each iteration of the method, the co-state initial value was quickly 
calculated with an efficient technique described in the Examples section. Fi- 
nally, the efficiency and reliability of the method were discussed, and the 
“dsolve” package of Maple 18 was used as the exact solution to check the 
validity of the method. 
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